%% Overview
% This file calculates all of the tables in the paper.
% Note that this file loads a pseudo data file instead of the true data
% file.


load pseudo_data_file;

%% Calculate persistences
% Earnings growth using 1yr and 3-5yr expectations
index_temp = e_growth_long(:,2)~=0;
[coeff,se_coeff] = beta_hac(e_growth_long(index_temp,2)*3,e_growth_long(index_temp,1));
temp = roots([1 1 1 0 -coeff]);
phi_e = real(temp(find(abs(imag(temp))==min(abs(imag(temp))),1)));
temp_high = roots([1 1 1 0 -coeff-se_coeff]);
phi_e_real_high = real(temp_high(find(abs(imag(temp_high))==min(abs(imag(temp_high))),1)));
temp_low = roots([1 1 1 0 -coeff+se_coeff]);
phi_e_real_low = real(temp_low(find(abs(imag(temp_low))==min(abs(imag(temp_low))),1)));
se_phi_e = abs(phi_e_real_high - phi_e_real_low)/2;


% Inflation using 1yr and 1-10yr expectations
index_temp = ~isnan(spf_final(:,2));
[coeff,se_coeff] = beta_hac(spf_final(index_temp,2)*10-spf_final(index_temp,1),spf_final(index_temp,1));
temp = roots([ones(9,1); -coeff]);
phi_pi = temp(imag(temp)==0);
grad = abs(sum((1:9).*(phi_pi.^(0:8))));
se_phi_pi = se_coeff/grad;


%% Figures 1 and 2
index_e_temp = e_growth_long(:,2)~=0;
figure; plot(sample_long,e_growth_long(:,1));
hold on; plot(sample_long(index_e_temp),e_growth_long(index_e_temp,2));
legend('Note that this is pseudo-data','location','best');

index_pi_temp = ~isnan(spf_final(:,2));
figure; plot(sample_long,spf_final(:,1));
hold on; plot(sample_long(index_pi_temp),spf_final(index_pi_temp,2));
legend('Note that this is pseudo-data','location','best');

%% Table 1: testing single condition
% This table tests the single necessary and sufficient condition for
% pe and the yield.

% 3-5 year earnings growth
e_growth_emp_long_3_5 = zeros(length(sample_long),1);
e_growth_emp_long_3_5(1:end-16) = (e_growth_emp_long(9:end-8) + e_growth_emp_long(13:end-4) + e_growth_emp_long(17:end))/3;
fe_e_3_5 = zeros(length(sample_long),1);
index_e_temp = e_growth_long(:,2)~=0 & e_growth_emp_long_3_5~=0;
fe_e_3_5(index_e_temp) = e_growth_emp_long_3_5(index_e_temp) - e_growth_long(index_e_temp,2);
index_pi_temp = ~isnan(spf_final(:,2)) & ~isnan(infl_emp_long_10);


table1 = zeros(6*3-1,3);
[table1(1,:), table1(2,:)] = cov_hac(100*[spf_final(:,1) infl_emp_long infl_emp_long-spf_final(:,1)],pe_long);
[table1(4,:), table1(5,:)] = cov_hac(100*10*[spf_final(index_pi_temp,2) infl_emp_long_10(index_pi_temp)  infl_emp_long_10(index_pi_temp)-spf_final(index_pi_temp,2)],pe_long(index_pi_temp));
[table1(7,:), table1(8,:)] = cov_hac(100*[spf_final(:,1) infl_emp_long infl_emp_long-spf_final(:,1)],yields_final(:,2));
[table1(10,:), table1(11,:)] = cov_hac(100*10*[spf_final(index_pi_temp,2) infl_emp_long_10(index_pi_temp)  infl_emp_long_10(index_pi_temp)-spf_final(index_pi_temp,2)],yields_final(index_pi_temp,2));
[table1(13,:), table1(14,:)] = cov_hac(100*[e_growth_long(:,1) e_growth_emp_long e_growth_emp_long-e_growth_long(:,1)],pe_long);
[table1(16,:), table1(17,:)] = cov_hac(100*3*[e_growth_long(index_e_temp,2) e_growth_emp_long_3_5(index_e_temp) fe_e_3_5(index_e_temp)],pe_long(index_e_temp));

%% Table 2: decomposing comovements
% Check that result for short-term earnings growth is mostly explained by
% the response to the surprise. Check that result for inflation is mostly
% explained by believed persistence of inflation.

eps_pi = spf_final(index_pi_temp,2)*10 - (1-phi_pi^10)/(1-phi_pi)*spf_final(index_pi_temp,1);
beta_e = 0.74;
eps_e = e_growth_long(5:end,1) + beta_e*(e_growth_emp_long(1:end-4) - e_growth_long(1:end-4,1));
table2 = zeros(8,3);
[table2(1,:), table2(2,:)] = cov_hac([10*spf_final(index_pi_temp,2)*100  (1-phi_pi^10)/(1-phi_pi)*spf_final(index_pi_temp,1)*100 eps_pi*100],pe_long(index_pi_temp));
[table2(4,:), table2(5,:)] = cov_hac([10*spf_final(index_pi_temp,2)*100  (1-phi_pi^10)/(1-phi_pi)*spf_final(index_pi_temp,1)*100 eps_pi*100],yields_final(index_pi_temp,2));
[table2(7,:), table2(8,:)] = cov_hac([100*e_growth_long(5:end,1) -beta_e*(e_growth_emp_long(1:end-4) - e_growth_long(1:end-4,1))*100 100*eps_e],pe_long(5:end));

%% Table 3: robustness to other price ratios

% Repeat the covariance tests using pd and pe_cape.
table3 = zeros(15*3-1,3);
table3(1:2,:) = table1(1:2,:);
[table3(4,:), table3(5,:)] = cov_hac(100*[spf_final(:,1) infl_emp_long infl_emp_long-spf_final(:,1)],pd_long);
[table3(7,:), table3(8,:)] = cov_hac(100*[spf_final(:,1) infl_emp_long infl_emp_long-spf_final(:,1)],pe_cape_long);
table3(10:11,:) = table1(4:5,:);
[table3(13,:), table3(14,:)] = cov_hac(100*10*[spf_final(index_pi_temp,2) infl_emp_long_10(index_pi_temp)  infl_emp_long_10(index_pi_temp)-spf_final(index_pi_temp,2)],pd_long(index_pi_temp));
[table3(16,:), table3(17,:)] = cov_hac(100*10*[spf_final(index_pi_temp,2) infl_emp_long_10(index_pi_temp)  infl_emp_long_10(index_pi_temp)-spf_final(index_pi_temp,2)],pe_cape_long(index_pi_temp));
table3(19:20,:) = table1(13:14,:);
[table3(22,:), table3(23,:)] = cov_hac(100*[e_growth_long(:,1) e_growth_emp_long e_growth_emp_long-e_growth_long(:,1)],pd_long);
[table3(25,:), table3(26,:)] = cov_hac(100*[e_growth_long(:,1) e_growth_emp_long e_growth_emp_long-e_growth_long(:,1)],pe_cape_long);

% Repeat the covariance test for short-term nominal earnings growth
% expectations using different denominators.
table3(28:29,:) = table1(13:14,:);
[table3(31,:), table3(32,:)] = cov_hac(100*[e_growth_long(:,1)+e_long-d_long e_growth_emp_long+e_long-d_long e_growth_emp_long-e_growth_long(:,1)],pd_long);
[table3(34,:), table3(35,:)] = cov_hac(100*[e_growth_long(:,1)+e_long-e_cape_long e_growth_emp_long+e_long-e_cape_long e_growth_emp_long-e_growth_long(:,1)],pe_cape_long);


% Repeat covariance test for long-term nominal earnings growth expectations
% using pd and pe_cape.
table3(37:38,:) = table1(end-1:end,:);
[table3(40,:), table3(41,:)] = cov_hac(100*3*[e_growth_long(index_e_temp,2) e_growth_emp_long_3_5(index_e_temp) fe_e_3_5(index_e_temp)],pd_long(index_e_temp));
[table3(43,:), table3(44,:)] = cov_hac(100*3*[e_growth_long(index_e_temp,2) e_growth_emp_long_3_5(index_e_temp) fe_e_3_5(index_e_temp)],pe_cape_long(index_e_temp));


%% Table 4: Asset prices
rho = exp(mean(pd_long)) / (1+exp(mean(pd_long)));
pe_mod = mean(pe_long) + (e_growth_long(:,1)-mean(e_growth_long(:,1)))/(1-rho*phi_e) ...
       - (spf_final(:,1)-mean(spf_final(:,1)))/(1-rho*phi_pi);
pd_mod = pe_mod + e_long - d_long;
pe_cape_mod = pe_mod + e_long - e_cape_long;
coeff_yield = (1-phi_pi^10)/(1-phi_pi)/10;
yields_mod(:,1) = spf_final(:,1)-mean(spf_final(:,1)) + mean(yields_final(:,1));
yields_mod(:,2) = coeff_yield*(spf_final(:,1) - mean(spf_final(:,1))) + mean(yields_final(:,2));
[pe_mod_lt,pe_mod_st] = hpfilter(pe_mod,1600); % HP filter
[pe_lt,pe_st] = hpfilter(pe_long,1600); % HP filter

% Table of regressions of asset prices
table4 = zeros(6*3-1,4);
[~,se,coeff] = hac(pe_mod,pe_long,'Display','off');
[~,~,r_square,r_square_res] = beta_gen(pe_long,pe_mod);
table4(1,:) = [coeff(1) coeff(2) r_square r_square_res];
table4(2,1:2) = se';

[~,se,coeff] = hac(pd_mod,pd_long,'Display','off');
[~,~,r_square,r_square_res] = beta_gen(pd_long,pd_mod);
table4(4,:) = [coeff(1) coeff(2) r_square r_square_res];
table4(5,1:2) = se';

[~,se,coeff] = hac(pe_cape_mod,pe_cape_long,'Display','off');
[~,~,r_square,r_square_res] = beta_gen(pe_cape_long,pe_cape_mod);
table4(7,:) = [coeff(1) coeff(2) r_square r_square_res];
table4(8,1:2) = se';

[~,se,coeff] = hac(yields_mod(:,2),yields_final(:,2),'Display','off');
[~,~,r_square,r_square_res] = beta_gen(yields_final(:,2), yields_mod(:,2));
table4(10,:) = [coeff(1) coeff(2) r_square r_square_res];
table4(11,1:2) = se';


[~,se,coeff] = hac(pe_mod_lt,pe_lt,'Display','off');
[~,~,r_square,r_square_res] = beta_gen(pe_lt,pe_mod_lt);
table4(13,:) = [coeff(1) coeff(2) r_square r_square_res];
table4(14,1:2) = se';

[~,se,coeff] = hac(pe_mod_st,pe_st,'Display','off');
[~,~,r_square,r_square_res] = beta_gen(pe_st,pe_mod_st);
table4(16,:) = [coeff(1) coeff(2) r_square r_square_res];
table4(17,1:2) = se';


% Test R^2 just using nominal earnings growth expectations
r_square_alt = zeros(3,1);
[~,~,~,r_square_alt(1)] = beta_gen(pe_long,e_growth_long(:,1)/(1-rho*phi_e));
[~,~,~,r_square_alt(2)] = beta_gen(pd_long,e_growth_long(:,1)/(1-rho*phi_e)+e_long-d_long);
[~,~,~,r_square_alt(3)] = beta_gen(pe_cape_long,e_growth_long(:,1)/(1-rho*phi_e)+e_long-e_cape_long);

%% Figure 3
figure; plot(sample_long, [pe_long pe_mod]); legend('Note that this is pseudo-data','location','best');

%% Table 5: Return predictability

% Calculate 10-year real and nominal returns.
r_nom_10 = zeros(length(sample_long),1);
for i = 1:10
    r_nom_10(1:end-36) = r_nom_10(1:end-36) + r_emp_long((i-1)*4+1:end-36+(i-1)*4);
end
r_real_10 = r_nom_10 - infl_emp_long_10*10;

table5 = zeros(8,7);
index = r_nom_10~=0;
[table5(1,1), table5(1,2), table5(1,3)] = beta_gen(r_nom_10(index),pe_long(index));
[table5(2,1), table5(2,2), table5(2,3)] = beta_gen(r_nom_10(index),pe_mod(index));
[table5(1,5), table5(1,6), table5(1,7)] = beta_gen(r_real_10(index),pe_long(index));
[table5(2,5), table5(2,6), table5(2,7)] = beta_gen(r_real_10(index),pe_mod(index));

[table5(4,1), table5(4,2), table5(4,3)] = beta_gen(r_nom_10(index),pd_long(index));
[table5(5,1), table5(5,2), table5(5,3)] = beta_gen(r_nom_10(index),pd_mod(index));
[table5(4,5), table5(4,6), table5(4,7)] = beta_gen(r_real_10(index),pd_long(index));
[table5(5,5), table5(5,6), table5(5,7)] = beta_gen(r_real_10(index),pd_mod(index));

[table5(7,1), table5(7,2), table5(7,3)] = beta_gen(r_nom_10(index),pe_cape_long(index));
[table5(8,1), table5(8,2), table5(8,3)] = beta_gen(r_nom_10(index),pe_cape_mod(index));
[table5(7,5), table5(7,6), table5(7,7)] = beta_gen(r_real_10(index),pe_cape_long(index));
[table5(8,5), table5(8,6), table5(8,7)] = beta_gen(r_real_10(index),pe_cape_mod(index));

%% Table 6: excess holding return volatility

% Calculate holding returns
r_emp_10yr = 10*yields_final(1:end-4,2) - 9*yields_final_9yr(5:end);
p_model_10 = -(1-phi_pi^10)/(1-phi_pi)*spf_final(:,1);
p_model_9 = -(1-phi_pi^9)/(1-phi_pi)*spf_final(:,1);
r_model_10yr = p_model_9(5:end) - p_model_10(1:end-4);

% Define variables for the variance ratios
% By defining new variables, we can detrend the variables or make other
% adjustments without writing over the original variables r_emp_10yr and
% r_model_10yr.
h_emp = detrend(r_emp_10yr);
r_emp = detrend(yields_final(:,1));
h_mod = detrend(r_model_10yr);
r_mod = detrend(yields_mod(:,1));


numlags = 12; % Number of lags to use if we assume variables are autocorrelated but not AR(1).
variables = {h_emp, r_emp, h_mod, r_mod};
variances = zeros(length(variables),1); % The variances for each variables accounting for the adjustment in effective observations
deg_freedom = zeros(length(variables),1); % The degrees of freedom for each variable.

% Check our assumption that autocorrelation is not significant after
% 12 lags. We do not find any significant autocorrelation.
b_check = zeros(length(variables),1);
se_check = zeros(length(variables),1);
for i = 1:length(variables)
    [b_check(i),se_check(i)] = persistence(variables{i},numlags);
end

% Calculate the variances, accounting for the adjusted number of
% observations.
for i = 1:length(variables)
    x = variables{i};
    T = length(x);

    autocorr_x = autocorr(x,'NumLags',numlags);
    coeff_aux = 1-(1:numlags)/T;
    T_eff = T/(1+2*(coeff_aux*autocorr_x(2:end)));
    deg_freedom(i) = T/(1+2*sum(autocorr_x(2:end).^2))-1;
    variances(i) = (T_eff/(T_eff-1))*var(x)*(T-1)/T;
end

% F-test of holding return variance
% We are testing the bound: Var(h) <= a*Var(r) where a = n^2/(2n-1)
n = 10;
a = n^2 / (2*n-1);

% Using CDF
p_val_emp = fcdf(variances(1)/(a*variances(2)),deg_freedom(1),deg_freedom(2),'upper');
p_val_mod = fcdf(variances(3)/(a*variances(4)),deg_freedom(3),deg_freedom(4),'upper');

table6 = [variances(1)/(a*variances(2)) variances(3)/(a*variances(4)); ...
          p_val_emp p_val_mod];
        

%% Fundamental extrapolation model
% Estimate objective values for persistences, variance and covariance of
% variables.
[phi_pi_obj, se_phi_pi_obj] = beta_hac(infl_emp_long,infl_emp_curr);
[phi_e_obj, se_phi_e_obj] = beta_hac(e_growth_emp_long,e_growth_emp_current_long);
[std_pi, se_std_pi] = std_est(infl_emp_long);
[std_e,se_std_e] = std_est(e_growth_emp_long);
[cov_pi_e,se_cov_pi_e] = cov_hac(e_growth_emp_long,infl_emp_long);

% Believed persistences already estimated earlier
phi_pi_star = phi_pi;
se_phi_pi_star = se_phi_pi;
phi_e_star = phi_e;
se_phi_e_star = se_phi_e;


% Estimate extrapolative weight (betas)
epsilon_pi = infl_emp_long-spf_final(:,1);
[beta_pi,se_beta_pi] = beta_hac(spf_final(5:end,1)-phi_pi_star*infl_emp_long(1:end-4),-epsilon_pi(1:end-4));
epsilon_e = e_growth_emp_long-e_growth_long(:,1);
[beta_e,se_beta_e] = beta_hac(e_growth_long(5:end,1)-phi_e_star*e_growth_emp_long(1:end-4),-epsilon_e(1:end-4));


% Table 7: model parameters
table7 = [phi_pi_obj phi_e_obj std_pi std_e cov_pi_e phi_pi_star phi_e_star beta_pi beta_e; ...
          se_phi_pi_obj se_phi_e_obj se_std_pi se_std_e se_cov_pi_e se_phi_pi_star se_phi_e_star se_beta_pi se_beta_e];




% Calculate diagnostic test in model
coeff_pi_y = (1-phi_pi_star^10)/(1-phi_pi_star)/10;
rho = exp(mean(pd_long)) / (1+exp(mean(pd_long)));
coeff_pi_pe = -1/(1-rho*phi_pi_star);
coeff_e_pe = 1/(1-rho*phi_e_star);

% Define variances and covariances in terms of variance of inflation and
% variance of nominal earnings growth.
c1_pi = (phi_pi_star-beta_pi)/(1-phi_pi_obj*beta_pi)* (phi_pi_star-beta_pi)*(1+phi_pi_obj*beta_pi)/(1-beta_pi^2); % Var(expected ST inflation)/Var(inflation)
c2_pi = (phi_pi_star-beta_pi)/(1-phi_pi_obj*beta_pi)*phi_pi_obj; % Cov(expected ST inflation, next year inflation) / Var(inflation)
c1_int = (phi_pi_star-beta_pi)*(phi_e_star-beta_e)/(1-beta_pi*beta_e)*(1/(1-phi_e_obj*beta_pi) + 1/(1-phi_pi_obj*beta_e) - 1); % Cov(expected ST inflation, expected ST nom e growth) / Cov(inflation, nom e growth)


var_pi_mod = std_pi^2;
var_e_mod = std_e^2;
cov_pi_e_mod = cov_hac(infl_emp_long,e_growth_emp_long);




c1_e = (phi_e_star-beta_e)/(1-phi_e_obj*beta_e)* (phi_e_star-beta_e)*(1+phi_e_obj*beta_e)/(1-beta_e^2); % Var(expected ST nom e growth)/Var(nom e growth)
c2_e = (phi_e_star-beta_e)/(1-phi_e_obj*beta_e)*phi_e_obj; % Cov(expected ST nom e growth, next year nom e growth) / Var(nom e growth)

diag_table = zeros(6,2);
diag_table(3,:) = [c1_pi c2_pi]*var_pi_mod*coeff_pi_y;
diag_table(4,1) = (1-phi_pi_star^10)/(1-phi_pi_star)*diag_table(3,1);
diag_table(4,2) = (1-phi_pi_obj^10)/(1-phi_pi_obj)*diag_table(3,2);
diag_table(1:2,1:2) = diag_table(3:4,1:2)*coeff_pi_pe / coeff_pi_y;
diag_table(5,:) = [(phi_e_star-beta_e)/(1-beta_e^2)*(1+phi_e_obj*beta_e) phi_e_obj] ...
                  *(phi_e_star-beta_e)/(1-phi_e_obj*beta_e)*var_e_mod*coeff_e_pe;
diag_table(6,1) = (phi_e_star^2+phi_e_star^3+phi_e_star^4)*diag_table(5,1);
diag_table(6,2) = (phi_e_obj^2+phi_e_obj^3+phi_e_obj^4)*diag_table(5,2);
diag_table = diag_table*100;

% Figure 4
figure;
subplot(3,1,1);
vals_data = [table1(1,1:2)'; table1(4,1:2)'];
vals_model = [diag_table(1,1:2)'; diag_table(2,1:2)'];
bar([vals_data vals_model]); 
ylim([-5 1]);
legend('Note that this is pseudo-data','location','best');
subplot(3,1,2);
vals_data = [table1(13,1:2)'; table1(16,1:2)'];
vals_model = [diag_table(5,1:2)'; diag_table(6,1:2)'];
bar([vals_data vals_model]);
ylim([0 15]);
subplot(3,1,3);
vals_data = [table1(7,1:2)'; table1(10,1:2)'];
vals_model = [diag_table(3,1:2)'; diag_table(4,1:2)'];
bar([vals_data vals_model]);
ylim([0 0.3]);

%% Replicating survey expectations (Figure 5)

eps_sim = zeros(length(sample_long),1);
eps_sim(1:4) = (phi_pi_star*infl_emp_curr(1:4) - spf_final(1:4,1))/beta_pi;
for i = 5:length(sample_long)
    eps_sim(i) = beta_pi*eps_sim(i-4) + infl_emp_curr(i) - phi_pi_star*infl_emp_curr(i-4); 
end
expect_sim = phi_pi_star*infl_emp_curr - beta_pi*eps_sim;

eps_sim_e = zeros(length(sample_long),1);
eps_sim_e(1:4) = (phi_e_star*e_growth_emp_current_long(1:4)-e_growth_long(1:4,1))/beta_e;
for i = 5:length(sample_long)
    eps_sim_e(i) = beta_e*eps_sim_e(i-4) + e_growth_emp_current_long(i)...
                 - phi_e_star*e_growth_emp_current_long(i-4); 
end
expect_sim_e = phi_e_star*e_growth_emp_current_long - beta_e*eps_sim_e;

eps_sim_e = zeros(length(sample_long),1);
temp = e_growth_long(:,1)-mean(e_growth_long(:,1));
eps_sim_e(1:4) = -temp(1:4)*2;
for i = 5:length(sample_long)
    eps_sim_e(i) = beta_e*eps_sim_e(i-4) + e_growth_emp_current_long(i)...
                 - phi_e_star*e_growth_emp_current_long(i-4); 
end
expect_sim_e = phi_e_star*e_growth_emp_current_long - beta_e*eps_sim_e;

figure; 
subplot(2,2,1); 
plot(sample_long,[expect_sim-mean(expect_sim) spf_final(:,1)-mean(spf_final(:,1))]);
subplot(2,2,2); 
index = ~isnan(spf_final(:,end));
coeff = (1-phi_pi_star^10)/(1-phi_pi_star)/10;
plot(sample_long(index),[coeff*(expect_sim(index)-mean(expect_sim(index))) spf_final(index,end)-mean(spf_final(index,end))]);
subplot(2,2,3);
plot(sample_long,[e_growth_long(:,1)-mean(e_growth_long(:,1)) expect_sim_e-mean(expect_sim_e)]);
ylim([-0.5 1.5]);
subplot(2,2,4);
index2 = e_growth_long(:,end)~=0;
coeff_e = (phi_e_star^2+phi_e_star^3+phi_e_star^4)/3;
plot(sample_long(index2),[coeff_e*(expect_sim_e(index2)-mean(expect_sim_e(index2))) e_growth_long(index2,end)-mean(e_growth_long(index2,end))]);
ylim([-0.5 1.5]);
legend('Note that this is pseudo-data','location','best');

%% Appendix: Alternative inflation expectations
% Compare our general results for SPF inflation expectations to the results
% using the Michigan inflation expectations. This includes the standard
% deviations, the subjective persistence, and the diagnostic test.
% Additionally, run the diagnostic test on the expectations from HPR.

% Compare standard deviations
index_scf_st = ~isnan(scf_final(:,1));
index_scf_lt = ~isnan(scf_final(:,2));
compare_std = [std(spf_final(:,1)) std(scf_final(index_scf_st,1)); std(spf_final(~isnan(spf_final(:,2)),2)) std(scf_final(index_scf_lt,2))];

% Compare the subjective persistence
% To calculate the subjective persistence, we treat the long-term Michigan
% expectations as 1-10yr expectations. In other words, we simply repeat the
% exact same process that we do for SPF but swap the short-term and
% long-term expectations with the Michigan data (without changing any of
% the horizons). 
[coeff_scf,se_coeff_scf] = beta_hac(scf_final(index_scf_lt,2)*10-scf_final(index_scf_lt,1),scf_final(index_scf_lt,1));
temp = roots([ones(9,1); -coeff_scf]);
phi_pi_scf = temp(imag(temp)==0);
grad_scf = abs(sum((1:9).*(phi_pi_scf.^(0:8))));
se_phi_pi_scf = se_coeff_scf/grad_scf;
compare_persistence = [phi_pi phi_pi_scf; se_phi_pi se_phi_pi_scf];


% Compare results of diagnostic test
index_scf_temp = ~isnan(scf_final(:,2)) & ~isnan(infl_emp_long_10);
table_michigan = zeros(11,3);
[table_michigan(1,:),table_michigan(2,:)] = cov_hac([infl_emp_long(index_scf_st) scf_final(index_scf_st,1) infl_emp_long(index_scf_st)-scf_final(index_scf_st,1)],pe_long(index_scf_st));
[table_michigan(4,:),table_michigan(5,:)] = cov_hac(10*[infl_emp_long_10(index_scf_temp) scf_final(index_scf_temp,2) infl_emp_long_10(index_scf_temp)-scf_final(index_scf_temp,2)],pe_long(index_scf_temp));
[table_michigan(7,:),table_michigan(8,:)] = cov_hac([infl_emp_long(index_scf_st) scf_final(index_scf_st,1) infl_emp_long(index_scf_st)-scf_final(index_scf_st,1)],yields_final(index_scf_st,2));
[table_michigan(10,:),table_michigan(11,:)] = cov_hac(10*[infl_emp_long_10(index_scf_temp) scf_final(index_scf_temp,2) infl_emp_long_10(index_scf_temp)-scf_final(index_scf_temp,2)],yields_final(index_scf_temp,2));
table_michigan = table_michigan(:,[2 1 3])*100; % Reorder columns and scale by 100 to match paper format for tables

%% Appendix: The role of LTG

e_emp_1_5 = [(e_growth_emp_long(1:end-16) + e_growth_emp_long(5:end-12) + e_growth_emp_long(9:end-8) ...
          +e_growth_emp_long(13:end-4)+e_growth_emp_long(17:end))/5; zeros(16,1)];
index = e_growth_long(:,end)~=0 & e_emp_1_5~=0;
ltg = e_growth_long(index,end);
realized = e_emp_1_5(index);

index1 = ismember(sample_long,bgls_mat(:,1));
index2 = ismember(bgls_mat(:,1),sample_long);
bgls = zeros(length(sample_long),1);
bgls(index1) = bgls_mat(index2,2);
bgls = bgls(index);

fe_bgls = realized - bgls;

% Test for overreaction in LTG
% Run test in levels (rather than logs) to exactly match the test in BGLS.
bgls_level = (exp(bgls)-1)*100;
fe_bgls_level = (exp(realized)-exp(bgls))*100;
[~,se,coeff] = hac([bgls_level(5:end)-bgls_level(1:end-4) bgls_level(1:end-4)],fe_bgls_level(5:end),'display','off');
overreaction_test_result = [coeff(2) se(2)];



% Matching the level of prices
p_long = pe_long + e_long;
kappa = log(1+exp(mean(pd_long))) - rho*mean(pd_long);
r_bar = mean(r_emp_long);
pe_perfect_foresight = zeros(length(p_long),1); % Calculate perfect foresight price
pe_perfect_foresight(end) = pe_long(end);
T = length(p_long);
rho_qua = rho^0.25;
for i = 1:T-1
    h = T-i;
    pe_perfect_foresight(i) = 0.25*(kappa-r_bar)*(1-rho_qua^h)/(1-rho_qua) + rho_qua^h*pe_long(end) + sum((e_long(i+1:end)-e_long(i:end-1)).*rho_qua.^(0:h-1)');
end
p_perfect_foresight = pe_perfect_foresight + e_long;


index_e = (e_growth_long(:,end)~=0);
p_syn_st = e_long(index_e) + e_growth_long(index_e,1);
temp = p_long(index_e) - p_syn_st;
[b,se] = beta_hac(temp,e_growth_long(index_e,end));
p_syn = e_long(index_e) + e_growth_long(index_e,1) + b*e_growth_long(index_e,end);
[b2,se2] = beta_hac(pe_long(index_e),e_growth_long(index_e,end));
p_syn3_lt = e_long(index_e) + b2*e_growth_long(index_e,end);
r_square_none = 1 - var(p_long(index_e) - e_long(index_e))/var(p_long(index_e));
r_square_st = 1 - var(p_long(index_e) - p_syn_st)/var(p_long(index_e));
r_square_full = 1 - var(p_long(index_e) - p_syn)/var(p_long(index_e));
r_square_lt = 1 - var(p_long(index_e) - p_syn3_lt)/var(p_long(index_e));
table_price_level = [1 1 1 1; 0 0 0 0; 1 1 0 0; 0 0 0 0; b 0 b2 0; se 0 se2 0; r_square_full r_square_st r_square_lt r_square_none];

% Standard deviation of price growth
p_long_alt = p_long(index_e);
p_perfect_foresight_alt = p_perfect_foresight(index_e);
table_price_growth_volatility = [std(p_long_alt(5:end)-p_long_alt(1:end-4)) std(p_perfect_foresight_alt(5:end)-p_perfect_foresight_alt(1:end-4)) std(p_syn(5:end)-p_syn(1:end-4)) std(p_syn_st(5:end)-p_syn_st(1:end-4))];



% Predicting returns with cash flow expectations
ed_exp = e_growth_long(:,1) + e_long - d_long;
e_growth_cape = e_growth_long(:,1) + e_long - e_cape_long;
exp_mat = [e_growth_long(:,1) ed_exp e_growth_cape e_growth_long(:,end)];
index = r_nom_10~=0;
index2 = index & e_growth_long(:,end)~=0;
table_ret_pred_appendix = zeros(4,4);
for i = 1:4
    if i == 4
        [b_nom,se_nom,r_square_nom] = beta_gen(r_nom_10(index2),exp_mat(index2,i));
        [~,~,~,r_square_res_nom] = beta_gen(r_nom_10(index2),-5*exp_mat(index2,i));
    else
        [b_nom,se_nom,r_square_nom] = beta_gen(r_nom_10(index),exp_mat(index,i));
        [~,~,~,r_square_res_nom] = beta_gen(r_nom_10(index),-1*exp_mat(index,i));
    end
    table_ret_pred_appendix(:,i) = [b_nom se_nom r_square_nom r_square_res_nom]';
end
table_ret_pred_appendix = table_ret_pred_appendix';




%% Appendix Robustness: value-weighted versus earnings-weighted

index1 = e_growth_long(:,end)~=0;
index2 = ismember(bgls_mat(:,1),sample_long(index1));
bgls_alt = bgls_mat(index2,1:2);
e_growth_long_alt = e_growth_long(index1,:);

% Volatility
table_value_earnings_weight = zeros(10,2);
table_value_earnings_weight(1,:) = [std(e_growth_long_alt(:,end)) std(bgls_alt(:,2))];

% Implied persistence
[coeff_bgls,se_coeff_bgls] = beta_hac(bgls_alt(:,2)*3,e_growth_long_alt(:,1));
temp = roots([1 1 1 0 -coeff_bgls]);
phi_e_bgls = real(temp(find(abs(imag(temp))==min(abs(imag(temp))),1)));
temp_high_bgls = roots([1 1 1 0 -coeff_bgls-se_coeff_bgls]);
phi_e_real_high_bgls = real(temp_high_bgls(find(abs(imag(temp_high_bgls))==min(abs(imag(temp_high_bgls))),1)));
temp_low_bgls = roots([1 1 1 0 -coeff_bgls+se_coeff_bgls]);
phi_e_real_low_bgls = real(temp_low_bgls(find(abs(imag(temp_low_bgls))==min(abs(imag(temp_low_bgls))),1)));
se_phi_e_bgls = abs(phi_e_real_high_bgls - phi_e_real_low_bgls)/2;
table_value_earnings_weight(3,:) = [phi_e phi_e_bgls];
table_value_earnings_weight(4,:) = [se_phi_e se_phi_e_bgls];


e_growth_emp_long_3_5 = zeros(length(sample_long),1);
e_growth_emp_long_3_5(1:end-16) = (e_growth_emp_long(9:end-8) + e_growth_emp_long(13:end-4) + e_growth_emp_long(17:end))/3;
e_growth_emp_long_3_5_bgls = e_growth_emp_long_3_5(index1);
fe_e_3_5_bgls = zeros(sum(index1),1);
index_e_temp_bgls = e_growth_emp_long_3_5_bgls~=0;
fe_e_3_5_bgls(index_e_temp_bgls) = e_growth_emp_long_3_5_bgls(index_e_temp_bgls)...
                                 - bgls_alt(index_e_temp_bgls,2);
pe_long_alt = pe_long(index1);
pd_long_alt = pd_long(index1);
pe_cape_long_alt = pe_cape_long(index1);
[table_value_earnings_weight(6:13,1)] = table3(37:44,1);
[table_value_earnings_weight(6,2), table_value_earnings_weight(7,2)] = cov_hac(100*3*bgls_alt(index_e_temp_bgls,2),pe_long_alt(index_e_temp_bgls));
[table_value_earnings_weight(9,2), table_value_earnings_weight(10,2)] = cov_hac(100*3*bgls_alt(index_e_temp_bgls,2),pd_long_alt(index_e_temp_bgls));
[table_value_earnings_weight(12,2), table_value_earnings_weight(13,2)] = cov_hac(100*3*bgls_alt(index_e_temp_bgls,2),pe_cape_long_alt(index_e_temp_bgls));

[table_value_earnings_weight(15:22,1)] = table3(37:44,3);
[table_value_earnings_weight(15,2), table_value_earnings_weight(16,2)] = cov_hac(100*3*fe_e_3_5_bgls(index_e_temp_bgls),pe_long_alt(index_e_temp_bgls));
[table_value_earnings_weight(18,2), table_value_earnings_weight(19,2)] = cov_hac(100*3*fe_e_3_5_bgls(index_e_temp_bgls),pd_long_alt(index_e_temp_bgls));
[table_value_earnings_weight(21,2), table_value_earnings_weight(22,2)] = cov_hac(100*3*fe_e_3_5_bgls(index_e_temp_bgls),pe_cape_long_alt(index_e_temp_bgls));


%% Appendix Robustness: rescaling test by variance
% This appendix table repeats the diagnostic test but reports regression
% betas instead of covariances.

table_norm = zeros(6*3-1,3);
[table_norm(1,:), table_norm(2,:)] = beta_hac([spf_final(:,1) infl_emp_long infl_emp_long-spf_final(:,1)],pe_long);
[table_norm(4,:), table_norm(5,:)] = beta_hac(10*[spf_final(index_pi_temp,2) infl_emp_long_10(index_pi_temp)  infl_emp_long_10(index_pi_temp)-spf_final(index_pi_temp,2)],pe_long(index_pi_temp));
[table_norm(7,:), table_norm(8,:)] = beta_hac([spf_final(:,1) infl_emp_long infl_emp_long-spf_final(:,1)],yields_final(:,2));
[table_norm(10,:), table_norm(11,:)] = beta_hac(10*[spf_final(index_pi_temp,2) infl_emp_long_10(index_pi_temp)  infl_emp_long_10(index_pi_temp)-spf_final(index_pi_temp,2)],yields_final(index_pi_temp,2));
[table_norm(13,:), table_norm(14,:)] = beta_hac([e_growth_long(:,1) e_growth_emp_long e_growth_emp_long-e_growth_long(:,1)],pe_long);
[table_norm(16,:), table_norm(17,:)] = beta_hac(3*[e_growth_long(index_e_temp,2) e_growth_emp_long_3_5(index_e_temp) fe_e_3_5(index_e_temp)],pe_long(index_e_temp));



%% Appendix: Generalized persistence structure

% Inflation using 1yr and 1-10yr expectations
index_temp = ~isnan(spf_final(:,end));
[coeff,se_coeff] = beta_hac(spf_final(index_temp,end)*10-spf_final(index_temp,1),spf_final(index_temp,1));
temp = roots([ones(9,1); -coeff]);
phi_pi = temp(imag(temp)==0);
grad = abs(sum((1:9).*(phi_pi.^(0:8))));
se_phi_pi = se_coeff/grad;


% Get b_21,b_22
index_temp = e_growth_long(:,end)~=0;
[~,se_temp,b_temp] = hac([spf_final(index_temp,1) e_growth_long(index_temp,1)],3*e_growth_long(index_temp,end),'display','off');
b_21 = b_temp(2);
b_22 = b_temp(3);

% Estimate phi_e_interact
temp = roots([1 1 1 0 -b_22]);
phi_e_interact = real(temp(find(abs(imag(temp))==min(abs(imag(temp))),1)));
temp_high = roots([1 1 1 0 -b_22-se_temp(3)]);
phi_e_real_high_interact = real(temp_high(find(abs(imag(temp_high))==min(abs(imag(temp_high))),1)));
temp_low = roots([1 1 1 0 -b_22+se_temp(3)]);
phi_e_real_low_interact = real(temp_low(find(abs(imag(temp_low))==min(abs(imag(temp_low))),1)));
se_phi_e_interact = abs(phi_e_real_high_interact - phi_e_real_low_interact)/2;

% Estimate phi_pi,e (interaction term)
coeff_temp = phi_pi+phi_pi^2+phi_pi^3 ...
           + phi_e_interact+phi_e_interact^2+phi_e_interact^3 ...
           + phi_pi*phi_e_interact+phi_pi^2*phi_e_interact+phi_pi*phi_e_interact^2;
phi_pi_e = b_21 / coeff_temp;
se_phi_pi_e = se_temp(2)/coeff_temp;

% Including the interaction term does not change the results for asset
% pricing or return predictability.
rho = exp(mean(pd_long)) / (1+exp(mean(pd_long)));
pe_mod_interact = mean(pe_long) + (e_growth_long(:,1)-mean(e_growth_long(:,1)))/(1-rho*phi_e_interact) ...
       - (spf_final(:,1)-mean(spf_final(:,1)))/(1-rho*phi_pi)*(1-rho*phi_pi_e/(1-rho*phi_e_interact));
[b,se,r_square,r_square_res] = beta_gen(pe_long,pe_mod_interact);
index = r_nom_10~=0;
[b_return, se_return, r_square_return] = beta_gen(r_real_10(index),pe_mod_interact(index));
temp = r_real_10 + infl_emp_long_10*10;
[b_return_nom, se_return_nom, r_square_return_nom] = beta_gen(temp(index),pe_mod_interact(index));

